#####Jill Shepperd and Jana von Stein. R code for "Attitudes and Action in International Refugee Policy:Evidence from Australia." 
#####Please direct questions to jana.vonstein@anu.edu.au 
library(foreign); library(survival); library(ggplot2); library(coefplot); library(ggpubr); library(mediation); library(MASS); library(survey); library(margins); library(effects); library(colorRamps)
library(readstata13); library(reshape2); library(lattice); library(scales); library(ggthemes); library(RColorBrewer); library(effects);library(allEffects); library(readstata13); library(stargazer); library(ATE)
library(expss);library(clubSandwich); library(EdSurvey);library(Rmisc);library(ggpubr);library("gridExtra"); library(jtools);library(allEffects)
library(factoextra); library(pwr); library(powerMediation);library(aod);library(WebPower)

setwd('~/<INSERT WORKING DIRECTORY HERE')
data<- read.dta13(file='shepperd-vonstein-final.dta', convert.underscore=TRUE)
##############FIGURES################
#############Figure 2
#if Figure 2 doesn't generate properly, reload the data using the code above. Sometimes R gets persnickety and turns everything to "NA"
#tell R treatment group is a factor
data$treatmentgroup <- factor(data$treatmentgroup, levels = c("control", "il", "moral", "reputational"),
  labels = c("Control", "Int'l Law","Moral", "Reputational"))     
#generate means and standard errors of policy approval for each treatment group
sum_policyapproval <- summarySE(data, measurevar =  "policyapproval1", groupvar = c("treatmentgroup"),
  na.rm = TRUE,conf.interval = .9)
#plot means and standard errors of policy approval by treatment group
ggplot(sum_policyapproval, aes(x = treatmentgroup, y = policyapproval1)) + 
  geom_bar(position=position_dodge(), stat="identity", fill = "grey50")+
  geom_errorbar(aes(ymin = policyapproval1 - ci, ymax = policyapproval1 + ci), width=.1) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank())+
  labs(x = "Treatment Group", y = "Mean Policy Approval") +
  coord_cartesian(ylim = c(2, 3))+theme( axis.line.y = element_line(colour = "grey50"))+
  theme( axis.line.x = element_line(colour = "grey50"))+ggtitle("Figure 2. Attitudes Toward Australia’s Boat Arrival Policy")+
  theme(plot.title = element_text(hjust = 0.5))
#################FIGURE 3
#OLS models of policy approval. These are just hypothesis tests, but in order to plot all three in one figure, we had to make replicate
#dummy variables for the treatment groups.
#intl law
policyapproval.ilcontrol.base<-lm(policyapproval1~ilframe, data=data,subset=ilframe==1|noframe==1);summary(policyapproval.ilcontrol.base)
policyapproval.ilmoral.base<-lm(policyapproval1~ilframe2, data=data,subset=ilframe==1|moralframe==1);summary(policyapproval.ilmoral.base)
policyapproval.ilreputational.base<-lm(policyapproval1~ilframe3, data=data,subset=ilframe==1|reputationalframe==1);summary(policyapproval.ilreputational.base)
policyapprovalil.plot<-plot_coefs(policyapproval.ilcontrol.base, policyapproval.ilmoral.base, policyapproval.ilreputational.base,
  scale = FALSE,inner_ci_level = .9, coefs=c("Control"="ilframe","Moral"="ilframe2", "Reputational"="ilframe3"),
  colors=c("dodgerblue2","dodgerblue2","dodgerblue2"),point.shape = FALSE)
policyapprovalil.plot<-policyapprovalil.plot+xlab("Average Treatment Effect")+ylab("") +theme_bw()+
  theme(legend.position="none")+ggtitle("Intl Law vs.")+theme(plot.title = element_text(hjust = 0.5))+
  theme(axis.text.y = element_text(size = 12))+theme(plot.margin=unit(c(.25,.1,.25,-.2), "cm"))+xlim(-.63,.48)
policyapprovalil.plot
#moral 
policyapprovalmoral.plot
policyapproval.moralcontrol.base<-lm(policyapproval1~moralframe, data=data,subset=moralframe==1|noframe==1)
policyapproval.moralil.base<-lm(policyapproval1~moralframe2, data=data,subset=moralframe==1|ilframe==1)
policyapproval.moralreputational.base<-lm(policyapproval1~moralframe3, data=data,subset=moralframe==1|reputationalframe==1)
policyapprovalmoral.plot<-plot_coefs(policyapproval.moralcontrol.base,policyapproval.moralil.base,policyapproval.moralreputational.base,
  scale = FALSE,inner_ci_level = .9, coefs=c("Control"="moralframe","Intl Law"="moralframe2", "Reputational"="moralframe3"),
  colors=c("dodgerblue2","dodgerblue2","dodgerblue2"),point.shape = FALSE)
policyapprovalmoral.plot<-policyapprovalmoral.plot+xlab("Average Treatment Effect")+ylab("") +theme_bw()+
  theme(legend.position="none")+ggtitle("Moral vs.")+theme(plot.title = element_text(hjust = 0.5))+
  theme(axis.text.y = element_text(size = 12))+theme(plot.margin=unit(c(.25,0,.25,-.2), "cm"))+xlim(-.63,.48)
policyapprovalmoral.plot
#reputational 
policyapproval.reputationalcontrol.base<-lm(policyapproval1~reputationalframe, data=data,subset=reputationalframe==1|noframe==1)
policyapproval.reputationalil.base<-lm(policyapproval1~reputationalframe2, data=data,subset=reputationalframe==1|ilframe==1)
policyapproval.reputationalmoral.base<-lm(policyapproval1~reputationalframe3, data=data,subset=reputationalframe==1|moralframe==1)
policyapprovalreputational.plot<-plot_coefs(policyapproval.reputationalcontrol.base,policyapproval.reputationalil.base,
  policyapproval.reputationalmoral.base, scale = FALSE,inner_ci_level = .9,
  coefs=c("Control"="reputationalframe","Intl Law"="reputationalframe2", "Moral"="reputationalframe3"),
  colors=c("dodgerblue2","dodgerblue2","dodgerblue2"),point.shape = FALSE)
policyapprovalreputational.plot<-policyapprovalreputational.plot+xlab("Average Treatment Effect")+ylab("") +theme_bw()+
  theme(legend.position="none")+ggtitle("Reputational vs.")+theme(plot.title = element_text(hjust = 0.5))+
  theme(axis.text.y = element_text(size = 12))+theme(plot.margin=unit(c(.25,.15,.25,-.2), "cm"))+xlim(-.63,.48)
policyapprovalreputational.plot
#generate final Figure 3
grid.arrange(policyapprovalil.plot,policyapprovalmoral.plot,policyapprovalreputational.plot,ncol=3)
#####################################################################################################################
setwd('~/Dropbox/australia-frames/final-submission')
data<- read.dta13(file='shepperd-vonstein-final.dta', convert.underscore=TRUE)
#############FIGURE 4
#if Figure 4 doesn't generate properly, reload the data using the code above. Sometimes R gets persnickety and turns everything to "NA"
#tell R treatment group is a factor
data$treatmentgroup <- factor(data$treatmentgroup, levels = c("control", "il", "moral", "reputational"),labels = c("Control", "Int'l Law","Moral", "Reputational"))     
#generate means and standard errors of policy approval for each treatment group
sum_petition <- summarySE(data, measurevar =  "petition1", groupvar = c("treatmentgroup"),na.rm = TRUE,conf.interval = .9)
#plot means and standard errors of petition by treatment group
petition.plot<-ggplot(sum_petition, aes(x = treatmentgroup, y = petition1)) + 
  geom_bar(position=position_dodge(), stat="identity", fill = "grey50")+
  geom_errorbar(aes(ymin = petition1 - ci, ymax = petition1 + ci), width=.1) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank())+
  labs(x = "Treatment Group", y = "% Responded Yes") +
  coord_cartesian(ylim = c(0, .4))+theme( axis.line.y = element_line(colour = "grey50"))+
  theme( axis.line.x = element_line(colour = "grey50"))+ggtitle("Sign Petition")+theme(plot.title = element_text(hjust = 0.5))
#plot means and standard errors of protest by treatment group
sum_protest<- summarySE(data, measurevar =  "protest1", groupvar = c("treatmentgroup"),na.rm = TRUE,conf.interval = .9)
#plot means and standard errors of policy approval by treatment group
protest.plot<-ggplot(sum_protest, aes(x = treatmentgroup, y = protest1)) + 
  geom_bar(position=position_dodge(), stat="identity", fill = "grey50")+
  geom_errorbar(aes(ymin = protest1 - ci, ymax = protest1 + ci), width=.1) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank())+
  labs(x = "Treatment Group", y = "% Responded Yes") +
  coord_cartesian(ylim = c(0, .4))+theme( axis.line.y = element_line(colour = "grey50"))+
  theme( axis.line.x = element_line(colour = "grey50"))+ggtitle("Protest")+theme(plot.title = element_text(hjust = 0.5))
#plot means and standard errors of protest by treatment group
sum_donate<- summarySE(data, measurevar =  "donate1", groupvar = c("treatmentgroup"),na.rm = TRUE,conf.interval = .9)
#plot means and standard errors of policy approval by treatment group
donate.plot<-ggplot(sum_donate, aes(x = treatmentgroup, y = donate1)) + 
  geom_bar(position=position_dodge(), stat="identity", fill = "grey50")+
  geom_errorbar(aes(ymin = donate1 - ci, ymax = donate1 + ci), width=.1) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank())+
  labs(x = "Treatment Group", y = "% Responded Yes") +
  coord_cartesian(ylim = c(0, .4))+theme( axis.line.y = element_line(colour = "grey50"))+
  theme( axis.line.x = element_line(colour = "grey50"))+ggtitle("Donate")+theme(plot.title = element_text(hjust = 0.5))
#generate Figure 4
grid.arrange(petition.plot,protest.plot,donate.plot,ncol=3)
##########FIGURE 5
#OLS models of policy action. These are just hypothesis tests, but in order to plot all three in one figure, we had to make replicate
#dummy variables for the treatment groups.

#intl law
factor.ilcontrol.base<-lm(factor~ilframe, data=data,subset=ilframe==1|noframe==1);summary(factor.ilcontrol.base)
factor.ilmoral.base<-lm(factor~ilframe2, data=data,subset=ilframe==1|moralframe==1);summary(factor.ilmoral.base)
factor.ilreputational.base<-lm(factor~ilframe3, data=data,subset=ilframe==1|reputationalframe==1);summary(factor.ilreputational.base)
factoril.plot<-plot_coefs(factor.ilcontrol.base, factor.ilmoral.base, factor.ilreputational.base,
  scale = FALSE,inner_ci_level = .9, coefs=c("Control"="ilframe","Moral"="ilframe2", "Reputational"="ilframe3"),
  colors=c("dodgerblue2","dodgerblue2","dodgerblue2"),point.shape = FALSE)
factoril.plot<-factoril.plot+xlab("Average Treatment Effect")+ylab("") +theme_bw()+
  theme(legend.position="none")+ggtitle("Intl Law vs.")+theme(plot.title = element_text(hjust = 0.5))+
  theme(axis.text.y = element_text(size = 12))+theme(plot.margin=unit(c(.25,.1,.25,-.2), "cm"))+xlim(-.22,.22)
factoril.plot

#moral
factor.moralcontrol.base<-lm(factor~moralframe, data=data,subset=moralframe==1|noframe==1)
factor.moralil.base<-lm(factor~moralframe2, data=data,subset=moralframe==1|ilframe==1)
factor.moralreputational.base<-lm(factor~moralframe3, data=data,subset=moralframe==1|reputationalframe==1)
factormoral.plot<-plot_coefs(factor.moralcontrol.base,factor.moralil.base,factor.moralreputational.base,
scale = FALSE,inner_ci_level = .9, coefs=c("Control"="moralframe","Intl Law"="moralframe2", "Reputational"="moralframe3"),
  colors=c("dodgerblue2","dodgerblue2","dodgerblue2"),point.shape = FALSE)
factormoral.plot<-factormoral.plot+xlab("Average Treatment Effect")+ylab("") +theme_bw()+
  theme(legend.position="none")+ggtitle("Moral vs.")+theme(plot.title = element_text(hjust = 0.5))+
  theme(axis.text.y = element_text(size = 12))+theme(plot.margin=unit(c(.25,.1,.25,-.2), "cm"))+xlim(-.22,.22)
factormoral.plot

#reputational
factor.reputationalcontrol.base<-lm(factor~reputationalframe, data=data,subset=reputationalframe==1|noframe==1)
factor.reputationalil.base<-lm(factor~reputationalframe2, data=data,subset=reputationalframe==1|ilframe==1)
factor.reputationalmoral.base<-lm(factor~reputationalframe3, data=data,subset=reputationalframe==1|moralframe==1)
factorreputational.plot<-plot_coefs(factor.reputationalcontrol.base,factor.reputationalil.base,
  factor.reputationalmoral.base, scale = FALSE,inner_ci_level = .9,
  coefs=c("Control"="reputationalframe","Intl Law"="reputationalframe2", "Moral"="reputationalframe3"),
 colors=c("dodgerblue2","dodgerblue2","dodgerblue2"),point.shape = FALSE)
factorreputational.plot<-factorreputational.plot+xlab("Average Treatment Effect")+ylab("") +theme_bw()+
  theme(legend.position="none")+ggtitle("Reputational vs.")+theme(plot.title = element_text(hjust = 0.5))+
  theme(axis.text.y = element_text(size = 12))+theme(plot.margin=unit(c(.25,.1,.25,-.2), "cm"))+xlim(-.22,.22)
factorreputational.plot
#Figure 5
grid.arrange(factoril.plot,factormoral.plot,factorreputational.plot,ncol=3)


#Table 1a, model 1
policyapproval.base<- lm(policyapproval1~ilframe+moralframe+reputationalframe, data=data); summary(policyapproval.base)
#Table 1a, model 2
policyapproval.extended<- lm(policyapproval1~ilframe+moralframe+reputationalframe+education1+income+agegroup+ethnicdist+
  female+labor+greens+otherindependent+ dontknowwontvote, data=data); summary(policyapproval.extended)
#output to Table 1a
stargazer(policyapproval.base,policyapproval.extended, ci=TRUE, type="html",out="action2.doc",star.cutoffs = c(.05, .01))

###Table 2a Wald tests###
linearHypothesis(policyapproval.base, "ilframe=moralframe")
linearHypothesis(policyapproval.base, "ilframe=reputationalframe")
linearHypothesis(policyapproval.base, "moralframe=reputationalframe")
linearHypothesis(policyapproval.extended, "ilframe=moralframe")
linearHypothesis(policyapproval.extended, "ilframe=reputationalframe")
linearHypothesis(policyapproval.extended, "moralframe=reputationalframe")

###Table 3a power tests
#subset treatment groups
control.subset<-subset(data,noframe==1)
il.subset<-subset(data,ilframe==1)
moral.subset<-subset(data,moralframe==1)
reputational.subset<-subset(data,reputationalframe==1)
#calculate N for each group
mean.controlpolicyapproval1<-mean(control.subset$policyapproval1)
n.control<-nrow(data[data$noframe == "1",])
n.il<-nrow(data[data$ilframe == "1",])
n.moral<-nrow(data[data$moralframe == "1",])
n.reputational<-nrow(data[data$reputationalframe == "1",])
#calculate mean for each group
mean.controlpolicyapproval1<-mean(control.subset$policyapproval1)
mean.ilpolicyapproval1<-mean(il.subset$policyapproval1)
mean.moralpolicyapproval1<-mean(moral.subset$policyapproval1)
mean.reputationalpolicyapproval1<-mean(reputational.subset$policyapproval1)
###Power analysis reported in Table 3a
pwr.t2n.test(n1 = n.il, n2=n.control , d = (mean.ilpolicyapproval1-mean.controlpolicyapproval1) , sig.level =.05 ); summary(n.il+n.control) #IL-control
pwr.t2n.test(n1 = n.il, n2=n.moral , d = (mean.ilpolicyapproval1-mean.moralpolicyapproval1) , sig.level =.05 );summary(n.il+n.moral)#IL-moral
pwr.t2n.test(n1 = n.il, n2=n.reputational , d = (mean.ilpolicyapproval1-mean.reputationalpolicyapproval1) , sig.level =.05 );summary(n.il+n.reputational)#IL-reputational
pwr.t2n.test(n1 = n.moral, n2=n.control , d = (mean.moralpolicyapproval1-mean.controlpolicyapproval1) , sig.level =.05 );summary(n.moral+n.control)#Moral-reputational
pwr.t2n.test(n1 = n.moral, n2=n.reputational , d = (mean.moralpolicyapproval1-mean.reputationalpolicyapproval1) , sig.level =.05 );summary(n.moral+n.reputational)#Moral-reputational
pwr.t2n.test(n1 = n.reputational, n2=n.control , d = (mean.reputationalpolicyapproval1-mean.controlpolicyapproval1) , sig.level =.05 );summary(n.reputational+n.control)#Reputational-control

###Table 5a, model 1
factor.base<- lm(factor~ilframe+moralframe+reputationalframe, data=data); summary(factor.base)
#Table 5a, model 2
factor.extended<- lm(factor~ilframe+moralframe+reputationalframe+education1+income+agegroup+ethnicdist+
  female+labor+greens+otherindependent+ dontknowwontvote, data=data); summary(policyapproval.extended)

stargazer(factor.base,factor.extended, ci=TRUE, type="html",out="action2.doc",star.cutoffs = c(.05, .01))

###Table 6a, Wald tests
linearHypothesis(factor.base, "ilframe=moralframe")
linearHypothesis(factor.base, "ilframe=reputationalframe")
linearHypothesis(factor.base, "moralframe=reputationalframe")
linearHypothesis(factor.extended, "ilframe=moralframe")
linearHypothesis(factor.extended, "ilframe=reputationalframe")
linearHypothesis(factor.extended, "moralframe=reputationalframe")

###Table 7a, power analysis
#calculate mean for each group
mean.controlfactor<-mean(control.subset$factor)
mean.ilfactor<-mean(il.subset$factor)
mean.moralfactor<-mean(moral.subset$factor)
mean.reputationalfactor<-mean(reputational.subset$factor)
###Power analysis reported in Table 3a
pwr.t2n.test(n1 = n.il, n2=n.control , d = (mean.ilfactor-mean.controlfactor) , sig.level =.05 ); summary(n.il+n.control) #IL-control
pwr.t2n.test(n1 = n.il, n2=n.moral , d = (mean.ilfactor-mean.moralfactor) , sig.level =.05 );summary(n.il+n.moral)#IL-moral
pwr.t2n.test(n1 = n.il, n2=n.reputational , d = (mean.ilfactor-mean.reputationalfactor) , sig.level =.05 );summary(n.il+n.reputational)#IL-reputational
pwr.t2n.test(n1 = n.moral, n2=n.control , d = (mean.moralfactor-mean.controlfactor) , sig.level =.05 );summary(n.moral+n.control)#Moral-reputational
pwr.t2n.test(n1 = n.moral, n2=n.reputational , d = (mean.moralfactor-mean.reputationalfactor) , sig.level =.05 );summary(n.moral+n.reputational)#Moral-reputational
pwr.t2n.test(n1 = n.reputational, n2=n.control , d = (mean.reputationalfactor-mean.controlfactor) , sig.level =.05 );summary(n.reputational+n.control)#Reputational-control

#########TABLES#################
###Table 8
#il group vs. control
ilcontrol.subset<-subset(data,ilframe==1|noframe==1)
model.m.ilcontrol<- lm(policyapproval1~ilframe, data = ilcontrol.subset)
model.y.ilcontrol<-lm(factor~ilframe+policyapproval1, data = ilcontrol.subset)
mediatedilcontrol.factor<- mediate(model.m.ilcontrol, model.y.ilcontrol, treat = "ilframe", mediator = "policyapproval1", boot=TRUE, sims=2000,conf.level = .95)
summary(mediatedilcontrol.factor);plot(mediatedilcontrol.factor)

#il group vs. moral
ilmoral.subset<-subset(data,ilframe==1|moralframe==1)
model.m.ilmoral<- lm(policyapproval1~ilframe, data = ilmoral.subset)
model.y.ilmoral<-lm(factor~ilframe+policyapproval1,  data = ilmoral.subset)
mediatedilmoral.factor<- mediate(model.m.ilmoral, model.y.ilmoral, treat = "ilframe", mediator = "policyapproval1", boot=TRUE, sims=1000)
summary(mediatedilmoral.factor); plot(mediatedilmoral.factor)

#il group vs. reputational
ilreputational.subset<-subset(data,ilframe==1|reputationalframe==1)
model.m.ilreputational<- lm(policyapproval1~ilframe, data = ilreputational.subset)
model.y.ilreputational<-lm(factor~ilframe+policyapproval1, data = ilreputational.subset)
mediatedilreputational.factor<- mediate(model.m.ilreputational, model.y.ilreputational, treat = "ilframe", mediator = "policyapproval1", boot=TRUE, sims=1000)
summary(mediatedilreputational.factor); plot(mediatedilreputational.factor)

#moral group vs. control
moralcontrol.subset<-subset(data,moralframe==1|noframe==1)
model.m.moralcontrol<- lm(policyapproval1~moralframe, data = moralcontrol.subset)
model.y.moralcontrol<-lm(factor~moralframe+policyapproval1,  data = moralcontrol.subset)
mediatedmoralcontrol.factor<- mediate(model.m.moralcontrol, model.y.moralcontrol, treat = "moralframe", mediator = "policyapproval1", boot=TRUE, sims=1000)
summary(mediatedmoralcontrol.factor);plot(mediatedmoralcontrol.factor)

#moral group vs. il group: no estimation -- this is just the inverse of il group vs. moral (above)

#moral group vs. reputational
moralreputational.subset<-subset(data,moralframe==1|reputationalframe==1)
model.m.moralreputational<- lm(policyapproval1~moralframe, data = moralreputational.subset)
model.y.moralreputational<-lm(factor~moralframe+policyapproval1,  data = moralreputational.subset)
mediatedmoralreputational.factor<- mediate(model.m.moralreputational, model.y.moralreputational, treat = "moralframe", mediator = "policyapproval1", boot=TRUE, sims=1000)
summary(mediatedmoralreputational.factor); summary(mediatedmoralreputational.factor)

#reputational group vs. control
reputationalcontrol.subset<-subset(data,reputationalframe==1|noframe==1)
model.m.reputationalcontrol<- lm(policyapproval1~reputationalframe, data = reputationalcontrol.subset)
model.y.reputationalcontrol<-lm(factor~reputationalframe+policyapproval1,  data = reputationalcontrol.subset)
mediatedreputationalcontrol.factor<- mediate(model.m.reputationalcontrol, model.y.reputationalcontrol, treat = "reputationalframe", mediator = "policyapproval1", boot=TRUE, sims=1000)
summary(mediatedreputationalcontrol.factor);plot(mediatedreputationalcontrol.factor)

#reputational group vs il group: no estimation -- this is just the inverse of il group vs. reputational (above)
#reputational group vs. moral group: no estimation -- this is just the inverse of moral group vs. reputational (above)

###Table 9a
#il-control
sigma.x.ilcontrol<-sd(ilcontrol.subset$ilframe)#calculate std dev of x for power analysis
sigma.m.ilcontrol<-sd(ilcontrol.subset$policyapproval1)#calculate std dev of mediator for power analysis
sigma.y.ilcontrol<-sd(ilcontrol.subset$factor)#calculate std dev of y for power analysis
summary(n.il+n.control)#sample size to report
a<- summary(model.m.ilcontrol)$coefficients[2] # path between predictor and mediator
b<- summary(model.y.ilcontrol)$coefficients[3] #path between mediator and dependent variable
wp.mediation(n = (n.il+n.control), power = NULL, a =a , b =b, varx =var(ilcontrol.subset$ilframe), vary=var(ilcontrol.subset$factor), 
  varm =var(ilcontrol.subset$policyapproval1), alpha = .05, interval = NULL)
               
#il-moral
sigma.x.ilmoral<-sd(ilmoral.subset$ilframe)#calculate std dev of x for power analysis
sigma.m.ilmoral<-sd(ilmoral.subset$policyapproval1)#calculate std dev of mediator for power analysis
sigma.y.ilmoral<-sd(ilmoral.subset$factor)#calculate std dev of y for power analyais
summary(n.il+n.moral)#sample size to report
a<- summary(model.m.ilmoral)$coefficients[2] # path between predictor and mediator
b<- summary(model.y.ilmoral)$coefficients[3] #path between mediator and dependent variable
wp.mediation(n = (n.il+n.moral), power = NULL, a =a , b =b, varx =var(ilmoral.subset$ilframe), vary=var(ilmoral.subset$factor), 
  varm =var(ilmoral.subset$policyapproval1), alpha = .05, interval = NULL)

#il-reputational
sigma.x.ilreputational<-sd(ilreputational.subset$ilframe)#calculate std dev of x for power analysis
sigma.m.ilreputational<-sd(ilreputational.subset$policyapproval1)#calculate std dev of mediator for power analysis
sigma.y.ilreputational<-sd(ilreputational.subset$factor)#calculate std dev of y for power analyais
summary(n.il+n.reputational)#sample size to report
a<- summary(model.m.ilreputational)$coefficients[2] # path between predictor and mediator
b<- summary(model.y.ilreputational)$coefficients[3] #path between mediator and dependent variable
wp.mediation(n = (n.il+n.reputational), power = NULL, a =a , b =b, varx =var(ilreputational.subset$ilframe), vary=var(ilreputational.subset$factor), 
  varm =var(ilreputational.subset$policyapproval1), alpha = .05, interval = NULL)

#moral-control
sigma.x.moralcontrol<-sd(moralcontrol.subset$moralframe)#calculate std dev of x for power analysis
sigma.m.moralcontrol<-sd(moralcontrol.subset$policyapproval1)#calculate std dev of mediator for power analysis
sigma.y.moralcontrol<-sd(moralcontrol.subset$factor)#calculate std dev of y for power analyais
summary(n.moral+n.control)#sample size to report
a<- summary(model.m.moralcontrol)$coefficients[2] # path between predictor and mediator
b<- summary(model.y.moralcontrol)$coefficients[3] #path between mediator and dependent variable
wp.mediation(n = (n.moral+n.control), power = NULL, a =a , b =b, varx =var(moralcontrol.subset$moralframe), vary=var(moralcontrol.subset$factor), 
  varm =var(moralcontrol.subset$policyapproval1), alpha = .05, interval = NULL)

###moral-reputational
sigma.x.moralreputational<-sd(moralreputational.subset$moralframe)#calculate std dev of x for power analysis
sigma.m.moralreputational<-sd(moralreputational.subset$policyapproval1)#calculate std dev of mediator for power analysis
sigma.y.moralreputational<-sd(moralreputational.subset$factor)#calculate std dev of y for power analyais
summary(n.moral+n.reputational)#sample size to report
a<- summary(model.m.moralreputational)$coefficients[2] # path between predictor and mediator
b<- summary(model.y.moralreputational)$coefficients[3] #path between mediator and dependent variable
wp.mediation(n = (n.moral+n.reputational), power = NULL, a =a , b =b, varx =var(moralreputational.subset$moralframe), vary=var(moralreputational.subset$factor), 
  varm =var(moralreputational.subset$policyapproval1), alpha = .05, interval = NULL)

##reputational-control
sigma.x.reputationalcontrol<-sd(reputationalcontrol.subset$reputationalframe)#calculate std dev of x for power analysis
sigma.m.reputationalcontrol<-sd(reputationalcontrol.subset$policyapproval1)#calculate std dev of mediator for power analysis
sigma.y.reputationalcontrol<-sd(reputationalcontrol.subset$factor)#calculate std dev of y for power analyais
summary(n.reputational+n.control)#sample size to report
a<- summary(model.m.reputationalcontrol)$coefficients[2] # path between predictor and mediator
b<- summary(model.y.reputationalcontrol)$coefficients[3] #path between mediator and dependent variable
wp.mediation(n = (n.reputational+n.control), power = NULL, a =a , b =b, varx =var(reputationalcontrol.subset$reputationalframe), vary=var(reputationalcontrol.subset$factor), 
  varm =var(reputationalcontrol.subset$policyapproval1), alpha = .05, interval = NULL)